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Abstract 

' We present the results of our lattice QCD study of the hadronic matrix elements 
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relevant to the physical radiative J ftp — > i] c j and h c — > r] c ^ decays. We used the 
twisted mass QCD action with iVf = 2 light dynamical quarks and from the com- 
putations made at four lattice spacings we were able to take the continuum limit. 
Besides the form factors parameterizing the above decays we also computed: (i) the 
hyperfine splitting and obtained A = 112 ± 4 MeV, (ii) the annihilation constant 



fj/ip which agrees with the one inferred from the measured T(J/tp — > e + e ). 
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1 Introduction 



After the observation of rjb at BaBar [T] the radiative decays of heavy quarkonia received a 
significant attention in the literature. As T(1S) — » 77^7 is not yet experimentally accessible 
due to the smallness of phase space, the experimenters turned to studying T(2S) —> r\bl 
and T(3S) — > 77^,7 modes from which they then extracted m Vb . A potential problem in 
that value is the insufficient control of theoretical uncertainties in the transition matrix 
elements to guarantee an accuracy at a percent level, especially when the radial excitations 
are involved. The corresponding transition matrix elements have been computed by using 
quark models [2]. Indeed the resulting value m^ b = 9390.9 ±2.1 MeV extracted from the 
BaBar experiment [1] was consistent with the value obtained from the similar measurements 
made at CLEO [I], leading to the following value of the hyperfine splitting [5], 

A£ xp - = m T(ls) - m Vb = 69.3 ± 2.8 MeV , (1) 

that turned out to be much larger than the values predicted by methods based on pertur- 
bative QCD, namely A b = 44 ± 11 MeV [6], 39 ± 14 MeV [7]. The lattice QCD results 
are inconclusive on this issue so far, although they seem to point towards the values larger 
than those obtained in refs. [Bl U\- For example, by using the simulations with Nf = 2 + 1 
staggered quarks and the Fermilab treatment of the heavy quarks on the lattice, the value 
Af u = 54 ± 12 MeV was obtained in ref. [8], while the non-relativistic QCD (NRQCD) 
treatment of the heavy quarks lead to A^ att = 70 ± 9 MeV [9]. Simulations of QCD with 
N{ = 2 + 1 light flavors of the domain wall light quark flavors and by using NRQCD for 
the heavy, resulted in Ajf tt = 60 ± 8 MeV [Hi]. 

The experimenters at Belle avoided using radial excitations and from a large sample of 
hb(lP) [TT] they were able to measure the hf,(lP) — > 772,7 decay rate, which resulted in a 
somewhat larger value m Vb = 9401. Oil. 9+H MeV pj] (i.e. smaller A° xp ' = 59.6±2.7 MeV), 
but still lighter than expected from the models [2] or the analytic calculations of refs j6j E]. 

Proponents of the extensions of the Standard Model involving more than one Higgs 
doublet speculated that the experimentally established pseudoscalar state 77^ might be 
actually a mixture of the true 77?, and the light parity-odd Higgs boson A pjj]. |^ This 
would solve the puzzle of too large a hyperfine splitting and would give more support to a 
plausible solution that m^o ~ 9 GeV. However, to give these speculations more support it 
is essential to check whether or not the hadronic matrix element used to extract m Vb from 
the mentioned experiments coincides with the results obtained by using the methods based 
on QCD from first principles. For example, by using NRQCD on the lattice, the authors of 
ref. [15] obtained much larger values for the transition matrix elements than those inferred 
from the measured T(nS) — > 77^7 (n = 2,3). Since the direct QCD simulations of the bb- 
systems are difficult because the lattice spacings are still too large to resolve the 6-quark 
mass, we decided to explore the similar physics processes in the charmed systems and 
study, J/ip — > r/ c 7 and h c (lP) — > r/ c 7. The established methodology of this paper will then 
be used for our future attempt to compute the amplitude for h b {\P) — > 77^7 decay on the 
lattice. Besides methodological issues these decays are physically interesting on their own. 

1 For a complete reviews with extensive lists of references please see ref. [H|3]. 

2 An abridged discussion on this issue with a more complete list of references can be found in ref. [Tl] . 
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One is the so called magnetic dipole (Ml) and the other electric dipole decay (El). Quark 
models fail to reproduce the measured Y(J/if) — > rj^) and, instead, obtain a significantly 
larger value [2]. On the other hand, h c (lP) has been discerned from the experimental 
background only recently and its dominant decay is indeed h c (lP) — > r) c "f, the branching 
fraction of which has been measured accurately. 

The first extensive study of the radiative decays of charmonia on the lattice has been 
reported in ref. [16] where the authors computed relevant matrix elements for a number of 
decay channels in the quenched approximation of QCD and with one lattice spacing. That 
computation has been extended to the case of iVf = 2 dynamical light quark flavors at 
single lattice spacing [17]. In this paper we will focus on J/if) — > r\ c ^ and h c (lP) — > -q^, for 
which we compute the desired form factors at four lattice spacings and then extrapolate 
them to the continuum limit. Those results may be used for a cleaner exclusion of the 
possibility of having very light m A o < 2m T (see e.g. [T3]). 



2 Hadronic Matrix Elements 

The transition matrix element responsible for the J /if) — > 7^7* decay reads, 

2 V(n 2 ) 

( Vc (k)\J™\J/i>(p, e A )> = eQ c W , (2) 

where J® m = Q c c7 M c is the relevant piece of the electromagnetic current, with Q c = 2/3 
in units of e = -y/47ra cm . Information regarding the non-perturbative QCD dynamics is 
encoded in the form factor V(q 2 ) and represents the most challenging part on the theory 
side. For the physical process, i.e. with the photon on-shell q 2 = 0, the decay rate is given 
by 



r(j/iMw) = % aeml J l \ 9 \v(o)\ 2 

27 (m J/v , + m Vi ^ 



2 

/ A \ 3 



8 / A \ 

= 7^ "em (mj/v, + m Vc ) \V(0) | z , (3) 

2 1 V m J/4' J 

where A stands for the hyperfine splitting A = mj/^ — m Vc . When both the initial and final 
hadrons are at rest the matrix element ([2]) is zero by definition. The smallest momentum 
that can be given to a hadron on the lattice with periodic boundary conditions is 2tt/L, 
which is very large for the lattices that we work with today and would make q 2 < 0, far 
from q 2 = 0. As a result we would have to work at several negative g 2 's, then model the q 2 
shape of the form factor as to extrapolate to the physical point, q 2 = 0. That methodology 
has been adopted in refs. [T6l [T7] . In this work, instead, we will use the so called twisted 
boundary conditions [18] which allow us to work directly at q 2 = 0. This is achieved by 
tuning the twisting angle 8q via the three momentum given to the pseudoscalar meson that 
fulfills the condition, 

m j/V _m ^ Lm^^-m 2 
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where we use q = (1, 1, 1) x 9 /L. For that purpose, and for each of our lattices, we first 
computed the masses of mj^ and of m^ c , and then by using eq. (jl]) we determined 9q that 
is then used in the computation of one of the charm quark propagators. In practice this 
last step is made by "twisting" the gauge links according to 

e*» /L UJx 



where 9^ = (0, 9) 



on which the quark propagator is computed according to, 

S°(x, 0; U) = e l ^ /L S c {x, 0; C/°) 



(5) 



(6) 



In our notation the quark propagator S c (x, 0) = 5 c (x, i; 0, 0) = (c(x)c(0))u, and we only 
in eq. ((6]) we write explicitly the gauge field configuration in the argument, to distinguish 
U from U e . In what follows U will be implicit. 

Similarly, in the computation of the physical h c — > rjcj decay, we compute the transition 
matrix element that is parameterized in terms of two form factors as, □ 



( Vc (k)\J™\h c (p,e x )) = ieQc { m^Ftf) I e 



The decay rate for the on-shell photon is [E] 
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(7) 



(8) 



To reach the physical form factor at q 2 = 0, with h c at rest, the twisted boundary condition 
applied on one of the charm quark propagators is made with 
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3 Two-point correlation functions 

Similarly to our recent publication [19], we use the gauge field configurations produced by 
ETM Collaboration [20J employing the maximally twisted mass QCD [2TJ, the details of 
which are summarized in tab. [TJ The masses of charmonia are extracted from the following 
correlation functions: 



(J2 Tr [S c (0, 0; x , t) l5 S' c (x, t; 0, 0) 75 ] ) , 

X 

1 3 



1=1 
3 



(10) 



3 From now on we will drop the label IP, and write h c only. 
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3.8 


3.9 


3.9 


4.05 


4.2 


4.2 


L 3 x T 


24 3 x 48 


24 3 x 48 


32 3 x 64 


32 3 x 64 


32 3 x 64 


48 3 x 96 


# meas. 


240 


240 


150 


150 


150 


100 


Mseal 


0.0080 


0.0040 


0.0030 


0.0030 


0.0065 


0.0020 


A*sea2 


0.0110 


0.0064 


0.0040 


0.0060 






Msea3 




0.0085 




0.0080 






Msea4 




0.0100 










a [fm] 


0.098(3) 


0.085(3) 


0.085(3) 


0.067(2) 


0.054(1) 


0.054(1) 


M<?o 2 ) P2] 


0.5816(2) 


0.6103(3) 


0.6103(3) 


0.6451(3) 


0.686(1) 


0.686(1) 


^(so 2 ) [22] 


0.746(11) 


0.746(6) 


0.746(6) 


0.772(6) 


0.780(6) 


0.780(6) 


123] 


0.2331(82) 


0.2150(75) 


0.2150(75) 


0.1849(65) 


0.1566(55) 


0.1566(55) 



Table 1: Summary of the details about the lattice ensembles used in this work (for more information 
see ref. [2D])- Data obtained at different /3's are rescaled by using the Sommer parameter ro/a, and the 
overall lattice spacing is fixed by matching f v obtained on the lattice with its physical value, leading to 
7'o = 0.440(12) fm (c.f. ref. [23.). All quark masses are given in lattice units. 



in which the Dirac structures are chosen to provide the coupling to the charmonium states 
with quantum numbers J PC = _+ , 1 , and l + ~, for rj c , J/ip and h c , respectively. 
S c (0,0; x,t) and S' c (0,0;x, t) refer to the propagators of the charm quark in the doublet 
ip{x) = [c(x) c'(x)] T entering the maximally twisted mass QCD action on the lattice [21] 



S = a 



X I fl 



9 /i /• 



+ // c U(z), (ii) 



and therefore the propagator S c (0,0; x,t) is obtained by inverting the above lattice Dirac 
operator with the Wilson parameter r, while S' c (0,0;x, t) is obtained by using — r. In 
practice r = 1 and m cr is the same as the one used in the production of the gauge field 
configurations [20]. Finally, V M and V* are the usual forward and backward derivatives 
on the lattice. In this study we also implement the Gaussian smearing on one of the 
sources [21]. In other words, one replaces c(x) by 



6k 



(12) 



4 Note that we write the action in the "physical basis" and not in the twisted one. 
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where the smearing operator H is defined via [25] 



■i 



(13) 



where is the n a times APE smeared link [26], defined in terms of (n a — 1) times smeared 



link U^ n ? ^ and its surrounding staples v} nc 



i, ii 



U na 



Projsu(3) 



wK-i) a (fto _i) 



(14) 



We chose the parameters 



4, n 5 



30, a = 0.5, n a = 20 



(15) 



which are kept fixed for all of our lattices. The value of the bare charm quark mass, fi c , at 
each of our lattices is given in tab. [1] It has been fixed according to the result of ref. [23] 
where it was shown that the charm quark computed from the comparison of the lattice 
results with the physical m Vc fully agrees with the value obtained by using the physical 
m Ds or m D . Therefore, we can say that m Vc , obtained by numerically solving 



(16) 



and then by fitting mf^(i) at large time separations to a constant, is merely a verification 
that, after a smooth continuum extrapolation, we indeed reproduce m^ p - = 2.980(1) GeV. 
To extract the values of mj/^ and rrif lc we proceed along the same line and compute 

m JMfcc(*) from 
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(17) 



In fig. [T] we show all three effective mass plots as obtained by using all four lattice spacings 
explored in this work and for one value of the sea quark mass which we choose to be the 
least light ones. We see that the effective masses for the pseudoscalar ad vector charmonia 
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Figure 1: Effective masses of the charmonium states, m® ff j,^ h (i), extracted from the two-point corre- 
lation functions according to eqs. (|16I17[) at four lattice spacings. Illustration is provided for one value of 
the sea quark mass. Note also that the smearing parameters used in this work are kept fixed to the same 
values for all our lattices. 



are excellent while the signal for the orbitally excited state, h c , is much more noisy. The 
effective masses are then combined to 



Rj/ip{t) 



mfit) 



RhM 



mf c (t) 



;is) 



The advantage of these ratios is that they have smaller statistical errors than any of the 
effective meson masses separately. In fig. [2] we show one such a ratio. On the plateaus, that 
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1.055 



^ 1.050 



1.045 




Figure 2: Plateau of Rj/^(t), defined in cq. (JTHJ) , obtained from our computations at (/3,/i sca ) = 
(4.2,0.0065). The shaded area is the result of the fit to a constant Rj/^. Results for all other values 
of (/3,/z S ea) explored in this study are listed in tab.[5J 



we carefully examined for each of our lattices, we then fit Rj/ip,h c (t) to a constant Rj/^^c- 
In tab. [2] we collect our results for Rj/^ t h c as obtained from all of the lattice ensembles at 
our disposal. As indicated in the plot in fig. [1] the fitting intervals involving the state h c 
are shorter. The shift of the plateau region to the right is made to account for the different 
lattice spacings, so that the fit is made at approximately the same physical separation 
between the interpolating field operators. 

After a smooth linear extrapolation to the continuum limit, 



R 



J/ip,h c 



pcont. 



1 + ^J/V,/ic m 9 + C J/^,h c 



(0.086 fm) : 



(19) 



we obtain 



= 1.0377(6) [ 1.0391(4) ] exp - , 
R c h ° nt - = 1.187(11) [ 1.1829(5) ] cxp - . (20) 

For the reader's convenience we also quoted the values obtained from experiments [5]. 
In eq. ( TT9"|) the parameter 6j/^,/ lc ~ measures the dependence on the sea quark mass, 
m q = m^ s (2 GeV), while the parameter cj/^^ c ~ 3 % measures the leading discretization 
effects. Division by 0/3=3.9 = 0.086 fm is made for convenience. The linear fit (Tl9l) describes 
our data very well except for the results obtained at (3 = 3.8. The results obtained at 
(3 = 3.8 can be either excluded from the continuum extrapolation, which is how we got 
the above results, or a term proportional to a 4 can be added in 019p which leads to a 
result that is fully consistent with the one quoted above. Finally, we should stress that the 
disconnected, OZI-suppressed, contributions to the correlation functions discussed in this 
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work have been neglected. The fact that our lattice results agree with the experimental 
values fl20|) can be viewed as a proof that the OZI suppressed contributions to the two-point 
functions are indeed very small. 



3.1 Dispersion Relation and the Sea Quark Mass Dependence 



We already mentioned that for the determination of the desired radiative decay form fac- 
tors one of the charm quark propagators is to be computed with twisted boundary condi- 
tions (jSJ), with the twisting angle tuned to ensure q 2 = 0, c.f. eqs. ( I4l9l) . In this section 
we check on the energy-momentum relation in the case of the pseudoscalar charmonium r\ c 
by exploring five different values of the twisting angle, covering the range of the meson's 
three-momenta, < \p\ < \J2 x 2it / L. We compute 



[E p ;t) = (^Tr[5f(0,0;f,t) 75 ^(f,t; 0,0) 75 ]) 



(21) 



with p — 9/L. For \8\ = we obviously get m Vc . For \6\ ^ 0, we proceed like in eq. ( fl6l) and 
fit Ep G (t) to a constant E p . As expected, due to the discretization effects the continuum 
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Figure 3: Test of the free boson lattice dispersion relation ([22]) on our results obtained at (3 = 4.05 in 
the case of the pseudoscalar charmonium j] c . 



relativistic formula, E 2 = m 2 c +p 2 , is not accurately verified on the lattice. Instead, our 
results satisfy 

4 sinh 2 ElL = 4 sin 2 M + C g 4 sinh 2 ^ , (22) 
2 2 1 2 y ' 

rather well. For Cp = 1 the above formula is the dispersion relation of a free boson on 
the lattice. After fitting our results to the above expression we find, for each of our lattice 



spacings, 



Cp = 0.97(1) 3 .8, 1.00(1)3.9, 1.03(1) 4 .05, 1.04(2) 4 .2 



(23) 



The results at /3 = 4.05 are shown in fig. [31 At this stage it is not clear whether or not 
with higher statistics all of the above Cr values would get closer to 1. The finite volume 



effects on the sea quark mass are unlikely to modify the expression (l22[) . That point we 
could check from our simulations at /3 = 3.9 and /i sea = 0.0040 where the results obtained 
at two lattices 24 3 x 48 and 32 3 x 64 are perfectly consistent. 

Since we are working with heavy quarks it is tempting to use the non-relativistic energy- 
momentum relation as well, but accounting for the lattice artifacts that according to ref. 



can be included by the distinction between the "rest" 
in, 



m 



) and the "kinetic" mass 



m 



Vc n (1) 



2m 



(24) 



</<- 



We fit our data to this expression and find that at each lattice spacing the kinetic mass 
is indeed larger than the rest one but they ultimately converge to the same value in the 
continuum limit, as they should on the basis of the restored Lorentz invariance, as shown 
in fig. HI Interestingly, however, we see that the discretization error in both the rest and 



h 1 r 



1 1 r 



"i 1 r 



3.6 
3.4 
3.2 
3.0 
2.8 
2.6 h 



m 



(i) 



4- 



m 



(0) 



_1 I L_ 



_1 I L_ 



_1 I L_ 







0.5 



(a/a 39 f 



(0) 



m^ c ' and m^J are the so called rest and kinetic mass obtained from the fit of our data to 



Figure 4: 

eq. (|24|) . The above plot shows extrapolation of their values (in physical units, GeV) to the continuum 
limit. 



kinetic masses are of the same size, but of the opposite sign. 



9 



0.20 



0.18 - 



0.16 



0.14 =~ 



0.12 



• 


(3=3.80 
p=3.90 


• 


▲ 


p=4.05 


< 


p=4-20 




0.06 



m q MS (2 GeV) 



Figure 5: Hyperfine splitting in charmonium as a function of the light sea quark mass [m s 
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MS 



(2 GeV)]. Four lines correspond to the fit A = A^ ' [l + 6a rn q ] made at each of our lattice spacings 



separately. All quantities are given in physical units [GeV] 



The above discussion on the dispersion relations is made on the data obtained at fixed 
value of the sea quark mass. We checked that indeed 

dm„. 



dm* 







(25) 



The same is, however, not true with mju, although this is hard to see from the mass ratio 
in eq. fTTTJT) alone. Instead we examined the hyperfine splitting, 



A = mj/^ - m Vc 
and from the fit of our data to, 



A = A 



cont. 



1 + b A m q + c A 



(0.086 fm) 5 



we find 



A 



cont. 



;il2±4)MeV [ 116.6 ± 1.2 ] cxp - 



(26) 



(27) 



(28) 



in good agreement with the experimental result written in brackets [5], and in excellent 
agreement with the result of lattice QCD computation presented in ref. [28], A = (111 ± 
5) MeV. Once more, this implicitly suggests that the contribution of the OZI breaking 
diagrams in the charmonia, neglected in our computations, are very small (c.f. discussion 
in ref. [29]). Note also that from the fit of our data to eq. (1271) we find, 



b A = 1.0(3) GeV" 1 , c A = 0.47(6) 



(29) 
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in qualitative agreement with ref. [8] where a tiny decrease of A is found while lowering 
the sea quark mass. Note, however, that this observation (pA ~ 0) disagrees with earlier 
findings of ref. [30]. To better appreciate that disagreement we fit our data at each lattice 
spacing separately as a linear function of the light sea quark mass and plot the resulting 
curves in fig. [5j We see that at each of our /3's the hyperfine splitting mildly decrease as the 
sea quark mass approaches the chiral limit. The above value for A cont ' is obtained without 
including the results at = 3.8 in the continuum extrapolation. If they are included our 
final result does not change but the error becomes smaller by 1 MeV. To be consistent with 
what we quote as a result of the continuum extrapolation in all other quantities computed 
in this work, we will quote A = (112 ± 4) MeV. 

4 Radiative Transition Form Factors 

To extract the desired hadronic matrix element (El) we computed the three point correlation 
functions 



= Tr [S' c (y; 0) 7 iS c (0, zfosZix, y)i 5 ] ) , (30) 

where P = cj^d ', Vi = crjid are the interpolating operators fixed at t = and t = t y = T/2 
(T being the time extension of our lattices). Using the fact that our three-momentum is 
isotropic, q = — (1, 1, 1) x 9 /L, we can average 

C v (q; t) = ~ [C 12 (q; t) + C 23 (q; t) + C 31 (q; t) - C 21 (q; t) - C 32 (q; t) - C 13 (q; t)} 
o 

^p. e -^ c (t flx -t) x go Mj/i, y(Q) x ^^ e _ mj/vjf ^ (31) 



where the last line is valid for the sufficiently separated operators in the correlation func- 
tion ( 1301) . which is ensured by the Gaussian smearing procedure. Coupling to the smeared 
pseudoscalar (vector) interpolating field operator P (V), is denoted by Zf, (Zy). [f| Elec- 
tromagnetic current, JJ m = Zylgfyc'jjC, is local and renormalized by using Zy{g^) listed 
in tab. [TJ Notice also that the pseudoscalar source operator has been fixed at t y = T/2, 
which simplifies the averaging of the signals propagating in both halves of the lattice, i.e. 
t G (0, ±T/2). In computing the propagators S c (x] y) we used the stochastic source tech- 
niques as explained in ref. [20]. Since the signals for the pseudoscalar and vector charmonia 
are very good, we can proceed to eliminate the sources in two ways. We can either divide 



5 The values of Zp and Zy are easily computed from the correlators of smeared source operators as 
C«(t) -> (\Z s P \ 2 /2m Vc ) exp[-ro,„f] , C^(t) -> (\Z V \ 2 /2m w ) exp[-m //v ,t] , 
to which we include the signal propagating from the opposite end of our lattice. 
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Figure 6: Plateaus exhibited by the semi- analytic and numerical ratios, i? S a,num(i) defined in eqs. (|32I33|) . 
are shown in the upper plots, while on the lower ones we show the semi-analytic ratios Rh(t) ()39[) . Since 
our results do not depend on the mass of the light sea quark, we plot the plateaus for the largest value of 
the sea quark mass at two values of /3. 



the correlator ( 130]) by the corresponding two point functions, 

C v (q;t) 



C^(t)CHE Q ,T/2-t) 



X ZpZy , 



(32) 



which we refer to as the numerical ratio, or by the analytic expression of the coupling of 
the source operators to the lowest states, 



Cy(q;t) 
Z S Z S 



(33) 



which we call the semi-analytic ratio. Obviously, the values for Zp V , mj/^ and E Vc are 
obtained from the study of the two-point correlation functions. The resulting plateaus 
for all the lattice spacings used in this work are shown in fig. [6j We see that thanks to 
efficiency of the smearing procedure the plateaus of the two ratios indeed coincide over 
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a broad range of time-slices. In the plots shown in fig. [6] we also indicate the intervals 
(shaded areas) in which we fitted the ratios -R sa ,num(^) to a constant. This is then used to 
obtain the form factor V(0), as indicated in eq. (j3~T]) . 

In tab. [2] we present our results for the form factor V(0), as obtained from all of the 
lattice data-sets and by using the ratio R sa .{t)- These results should be now extrapolated 
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Table 2: Detailed results for the hadronic quantities discussed in this paper, computed on each lattice 
data set specified in tab.Q] Note that m, (c and fj/j, are given in lattice units. 



to the physical limit (m sca = m q — )■ 0, a — )■ 0) by 



V(0) = V(0) 



cont. 



1 + b v m q + cy 



(0.086 fm) 5 



(34) 



from which we then obtain the two following values: 

1.937(34) (with (3 = 3.8) , 

1.941(35) (without (3 = 3.8) 



V(0) 



(35) 
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depending on whether or not we include the results obtained on our coarsest lattices in the 
continuum extrapolation. A few more comments are in order: 

• The values of the J/tp — > rjcj form factor computed on our lattices do not depend on 
the light sea quark mass, by ~ 0. 

• Although the smallness of discretization errors in tmQCD is guaranteed by construc- 
tion [they are 0(a 2 )}, they are quite large for the form factor V(0). More specifically, 
from the fit of our data to eq. fl34"|) we get c v = —23(2)%. Note however that in 
this case the results obtained at = 3.8 do not modify the result of the continuum 
extrapolation, i.e. the formula ( 134]) adequately describes all our results for V(0). 
This can be appreciated from the plot provided in fig. [TJ 

• The above results are obtained by using the semi-analytic ratio R sa . By repeating 
the same analysis with -R n um we verify the same features concerning the 0(a 2 ) effects 
and after extrapolating the results computed on all our lattices we obtain V(0) = 
1.903(48), entirely consistent with the results quoted in eq. ( 135]) . 

Our final result is: 

V(0) = 1.94(3) (1°) . (36) 

where the second error is a difference between the central values obtained by using the 
semi-analytic and numerical ratios discussed above. 
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Figure 7: Continuum extrapolation of the form factors V(0) and -Fi(O) computed on our lattices at 4 
lattice spacings. Note that for every value of a 2 we have several data points obtained for different value 
of the light sea quark mass. The continuum extrapolation is made according to eq. and similarly for 
the form factor -Fi(O). 
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We now turn to the discussion of the form factor Fi(0), relevant to the h c — > rjcj decay, 
as defined in eq. ([7]). To that end we compute the three point correlation functions 



C ijk (q;t) = X;<4(0) JT(x)P(v)) j qJ3 - yl 

= -(£ Tr [S' c (y; 0) Ti7 A(0, x) lk Sl{x, y) l5 ] ) . (37) 

Keeping in mind that we consider h c at rest and for q = —(1, 1, 1) x 6q/L, with 9q tuned 
as indicated in eq. fl9]), we can easily isolate the term proportional to -Fi(O) by combining: 

CM t) = ~ [C 123 (q; t) + C 231 (q; t) + C 312 (q; t)} , 

C (l; t) = \ [Cisi(g; t) + C 212 (q- 1) + C 323 (g; t) 

o 

+C 232 {q; t) + C 313 (q; t) + C 121 (q; t)} , 
C Fl (q;t) = C d (q;t) - C (q;t) 

z s , , , z s 



2E Vc " c ±v ' 2m t 



where, as before, in the last line we show the result of the spectral decomposition when all 
operators are sufficiently separated. By fitting 

*h(f) = lm[C jfs t)] >< ^ h E Vc e^WM , (39) 

to a constant we obtain m/ lc F 1 (0). The results for the form factor -Fi(O), as obtained from 
all of our lattice data sets, are presented in tab. [2j These results too need to be extrapolated 
to the continuum and chiral limits, in a way similar to eq. ( |34|) . 

• The error on -Fi(O) computed on each of our lattices is under 10%, but is nevertheless 
twice larger than those we have when computing ^(0). This is expected as the 
signal for h c is much harder to tame. For the same reason we could not compute the 
form factor by employing the numerical method, i.e. by dividing by the two-point 
correlation functions. 

• Contrary to the case of V(0), the discretization effects on -F\(0) are small (c.f. fig. E]). 
From the fit of our data to the form similar to eq. (1341) we find Cp 1 ~ 2%. We get: 

-0.57(2) (with = 3.8) , 

(40) 

-0.57(3) (without (3 = 3.8) , 
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Like in the case of V(0), within the accuracy of our data, the form factor Fi(0) is 
insensitive to the variation of the light sea quark mass. 



• Our final result is 

Fi(0) = -0.57(2) (1), (41) 

where the second error is our estimate of the uncertainty due to the method for 
extracting the form factor from the correlation functions. In the case of V(0) we were 
able to estimate that error from the difference between the central values obtained 
by using the semi-analytical ( 1331) and numerical ratios ( 1321) . That difference turned 
out to be less than 2%. We add the same error to i*i(0) leaving its sign free. 



4.1 J/ifj annihilation constant 

In addition to the above quantities we also computed the annihilation constant fj/^ defined 



as, 



(0|c(0)7„c(0)| JA% A)) = f m m m el , (42) 

which enters decisively in the expression for the well measured electronic width T(J/tp — > 
e + e~). Above, stands for the polarization vector of J/ip. fj/$ is computed along the 
same lines discussed in our previous paper [19], namely from the fit of the correlation 
function C^(t) from eq. (flOl) to the form 

cS*(t) m C ° S " |mw(r/2 '' )l Z^)(0|c(0) 7iC (0)|JM0,A)} V- , (43) 

where the non-perturbatively determined Za^q) are those listed in tab. [TJ In practice 
we of course combine the correlation functions in which both interpolating operators are 
smeared and those in which one operator is local and the other one is smeared. The results 
are given in tab. [2J which after extrapolating to the continuum by using a form analogous 
to eq. (jMD lead to 

/ w =414±8lgMeV. (44) 

The continuum extrapolation is smooth (c/ = 5(2)%) and the result is obtained by using 
all our lattices. If we leave out the results obtained at (3 = 3.8, the resulting fj/^ is larger 
by 9 MeV, which is the second error quoted above. Note also that fj/^ depends very mildly 
on the light sea quark mass (bf = 0.5(2) GeV -1 ). 



6 Note that in the three point functions we use t/* c wn i cn i s invariant under the same axial rota- 

tions that leave the twisted mass QCD action invariant and is the invariant vector current multiplicatively 
rcnormalized by Zv(go)- In the two point function, instead, Tpc^Ji^c (a = 1,2) is used which, at the 
maximal twist, corresponds to the axial current (in the twisted-mass basis) and is therefore rcnormalized 
by Z A (gl) m- 
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5 Phenomenology 



5.1 Decays of J/ip 

By using our result (144]) we can compute the electronic width of J/ip as 





cxp. 



(45) 



where we used a em = 1/134 [31], symmetrized the error bars, and quoted the experimentally 
established result [5] to better appreciate the agreement between our lattice result and 
experiment. Concerning the radiative decay J/ip — > rjcj, by inserting our value fl36|) in 
eq. OH]) we get 



where for the physical result we used the measured Br(J/ , — > ^7) = (1.7 ± 0.4)% 
and the full width Tj/^ = 92.9 ± 2.8 keV [5], as well as the physical values of m;^ = 
3096.92(1) MeV, and A = 116.6 ± 1.2 MeV. Had we used our A = 112 ± 3 MeV, the 
resulting T(J/ip — > ^7) would have been 10% smaller. 

Although somewhat lower than the quark model result, T(J/tp — > r/ c 7) = 2.85 keV j2], 
our lattice result obviously gives a larger J/ip — > rjcj decay rate, and the agreement with 
experiment is only at 2a. The effective theory approach of ref. [32] and the QCD sum rule 
analyses [33] succeeded at getting lower value for the decay rate of this decay, but with large 
uncertainties. Note that the dispersive (model independent) approach of ref. [31] predicted, 
T(J/ip — > rjc'j) = 2.2 -i- 3.2 keV, many years ago. All these results agree with ours too, 
except that we have smaller and controlled uncertainties. We hope more effort on the 
experimental side will be devoted to clarify the disagreement among various experiments, 
including the recent ones. For example, a study of this decay at BESIII would give us 
a very valuable information. Recent result at KEDR suggested a larger value for the 
branching fraction Br(J '/if) — > 7^7) = (2.34 ± 0.15 ± 0.40) % [35], which would result in 
T(J/i() — >• 7^7) = (2.2 ± 0.6) keV, in very good agreement with our result ( 1461) . 

As for the other lattice calculations of this form factor, we note that the quenched 
result of ref. [16], V(0) = 1.85(4), is only slightly lower than ours, while the one obtained 
at single lattice spacing with Nf = 2 light flavors in ref. [IT], 1^(0) = 2.01(2), is larger 
than our values at (3 = 4.05 listed in tab. [2j Apart from different methodology, a notable 
difference is that the authors of ref. [T7] used the point-split electromagnetic current that 
does not require renormalization, whereas we use the local current on the lattice that is 
properly renormalized by non-perturbatively determined Zylg^). Keep in mind that our 
final results is obtained after the extrapolation to the physical limit of the results obtained 
at several lattice spacings and for several different values of the light sea quark mass. 

5.2 h c — > rjc'j 

h c escaped the experimental detection for a long time and only recently CLEO succeeded 
to isolate this state [36] and observed that its prominent mode is precisely h c — > rjcj, the 



r(J/i[) r] c j) = 2.64(H)(3) keV 



[1.58(37) keV] 



exp. 



(46) 
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branching fraction of which was later accurately measured at the BESIII experiment, with 
a result: Br(h c — > r) c y) = (53 ± 7)% [3?]. We obviously cannot compute the branching 
ratio on the lattice, but with our form factor result ( 14TI) we can compute the decay width 
using eq. (JH]). We get 

T(h c ncj) = 0.72(5) (2) MeV . (47) 

This result can be combined with the measured Br(/i c — > n^) to estimate the width of the 
h c state. We obtain: 

r fe = — Vcl \ = 1.37 ±0.11 ±0.18 MeV, (48) 

Br(h c -)■ ?7 C 7) 

where the first error comes from our determination of the form factor -Fi(O), and the 
second one reflects the experimental uncertainty in the branching ratio. Notice also that 
we symmetrized the error bars. This constitutes a prediction that would be interestingto 
check against the actual experimental measurement once the latter becomes available. LJ 

To compare our estimate -Fi(O) = —0.57(2) with other lattice results we convert the 
value of the value reported in ref . [H] to our dimensionless form factor and obtain Fx (0) = 
—0.53(3), which agrees very well with our result. |f| Similar conversion of the result of 
ref. [T7] would result in -Fi(O) = —0.33(1), much smaller value than ours, whether we 
compare it with the values we obtain at f3 = 4.05 or the one in the continuum limit. 



6 Summary 

In this paper we presented results of our analysis of the radiative decays of charmonia 
by means of QCD simulations on the lattice. Using the twisted mass QCD with Nf = 2 
dynamical flavors at several small lattice spacings we were able to smoothly extrapolate 
the relevant form factors to the continuum limit. 

• We used the twisted boundary conditions to make sure that we extract the physical 
form factors, i.e. at q 2 = 0. We checked that at every lattice spacing explored in this 
paper our data indeed reproduce the latticized energy-momentum relation given in 
eq. (|2"2"|) . We also showed that the so-called kinetic and rest masses converge to the 
same value in the continuum limit, but that at fixed lattice spacing they both have 
discretization errors that are equal in size but different in sign. 

• We computed the hyperfine splitting and obtained 

A = m J/f -m Vc = 112 ±4 MeV, (49) 

and showed that it depends very mildly on the sea quark mass, with a slope being 
positive. From our computation of Rh c in eq. ( 120]) we obtain rrih c = 3.537(32) GeV, 
in good agreement with rn^' = 3.525 GeV. 

7 It will also be interesting to see the physics results of the effective theory developed in ref. [55] , 
8 More specifically, the relation between our i*i(0) and -E'i(O), defined in ref. [TT5] is: -Fi(O) = 
a t Ei (0)/m hc , with a t = 6.05(1) GeV, and a t Ei(0) = -0.306(14). 
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• We computed the hadronic form factor relevant to the J/x/j — > rjcj (Ml) transition 
and found 

V(0) = 1.92(3) (2), (50) 

which is larger than the one we would infer from the measured T(J/ip — > r) c ^f), 
although compatible at the 2a level. We found that the discretization effects are 
large and negative, but that they are adequately described by the linear function in 
a 2 . 

• Our result for the h c — > r]^ (El) transition form factor is 

F 1 (0) = -0.57(2)(l), (51) 

which mildly depends on the lattice spacing. Within the uncertainties quoted above, 
both our form factors are insensitive to the change of the sea quark mass. After 
combining our -Fi(O) with the measured Br(/i c — > rjc'-f) we deduced the value of the 
width, T hc = 1.37(22) MeV. 

• In addition to the above, we also computed the annihilation constant 

= 418 ± 8 ± 5 MeV , (52) 

which agrees with the measured decay width T(J/ip — > e + e~). 

In the above results we do not make any estimate of the size of systematic uncertainty due to 
the omitted s and c quarks in the sea. In ref. [39] it was claimed that the contributions from 
dynamical charm might be important. That point will be numerically assessed from the 
analysis similar to the one presented in this paper but on the set of gauge field configurations 
that include TVf = 2 + 1 + 1 dynamical quark flavors. We also emphasize that our results 
are obtained without inclusion of the OZI suppressed contributions. Their impact appears 
to be small in J/ip — > e + e~ decay, but their size in the radiative decays is unknown. 
They were neglected and the associated uncertainty cannot be estimated without actually 
attempting to compute the corresponding disconnected diagrams on the lattice. Such a 
computation would be very welcome. 

The strategy employed in this paper can be applied to compute the much needed 
F^ b ^ r,b (0) which we plan to do in the near future. 
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